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NUMERICAL ANALYSIS OF A PAIR OF SINGULAR INTEGRAL 
EQUATIONS WITH A PRINCIPAL VALUE 
by Raymond C. Cierkin 
Lewis Research Center 

SUMMARY 

This report describes a numerical analysis of a pair of integral equations. One 
equation is of the Volterra type and has a singularity which requires special handling. 

The other equation is of the Fredholm type with a principal value . 

The technique of integration of the principal value integral is described and tested on 
an improper integral with a known solution. To attain convergence of the iterative solu- 
tion by the method of Picard it was necessary to take two actions to overcome instability 
near the singularity: (1) reduce the relaxation parameter, and (2) constrain one function 
to behave smoothly near the singularity. 


INTRODUCTION 

Research on the shape of a magnetically balanced arc reported in reference 1 re- 
quired the solution of a pair of integral equations. Work done on the numerical solution 
of this problem has two features of interest to numerical analysts with similar problems: 
(1) the technique of handling the principal value of an integral, and (2) the method of 
achieving a satisfactory convergence of the iterative solution of the simultaneous integral 
equations despite the presence of a singularity. 

Because numerical analysts occasionally have need to solve problems involving prin- 
cipal values of integrals or to deal with integral equations having singular points, the ex- 
perience described in this report should prove helpful. 



SYMBOLS 


a 2 ,b 2 ,c 2 


D(s) 

f l’ f 2 

g(s) 


*3 

J(s) 

j 

^max 

K 

N 

P.V. 


q 

s. 

s 



X 


a v a 2 

P 

e 

e(s) 

e K (s) 

e 0 (s) 

r(s) 


coefficients in definition of f^, eq. (20) 
coefficients in definition of fg, eq. (21) 
difference in old and new 6 functions, eq. (5b) 
quadratic approximations used in eq. (19) 

known function of s used as test of accuracy of numerical method used to 
solve for 0(s) 

integral to left of principal value integral 

principal value integral 

integral to right of principal value integral 

approximating function to numerator integral in eq. (3) 

index number 

maximum index for 

refers to the iteration number 

number of points in interval 0^s<1.0 such that As = 1/(N - 1) 
principal value 

dummy variable of integration 

integral sum symmetric about q = s (eq. (16)) 

parameter 

t function j steps to left of s = 1 

11. 

function 9 on K m iteration if there were no relaxation 
dummy variable of integration in eq. (19) 
limits of integration in eq. (19) 
relaxation parameter 

variable approaching zero in defining 6(s) as principal value of integral 
one of two unknown functions of s 

1U 

function 9{ s) at K tn iteration 
initial guess at function 6(s) 
one of two unknown functions of s 
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r* (q) function r obtained by expanding r-^-(q) in first three terms of its Taylor 

series about q = s 

r’(s) derivative of r with respect to s 

t"(s) second derivative of r with respect to s 

r (q) function T(q) obtained by replacing t* by its finite difference analog 

a. 

th 

r„(s) function r(s) on K n iteration 

K. 

it. 

Tg(s) derivative of r with respect to s on K L iteration 

th 

Tj^( s) second derivative of r with respect to s on K n iteration 


STATEMENT OF THE PROBLEM 

The equations to be solved are 



where P.V. is the principal value, and the right side of equation (2) is defined as 




fs-e 

f 1 

lim 

1 1 

I r(q)dq f 1 j 

I T(q)dq 

€ -■*' 0 

7T 1 

1 q - s it 1 

1 q - s 


J 

-1 J 

s+e 


The quantities t(s) and 6( s) are the unknown functions over the interval 0 < s < 
t(-s) = r(s) defines t in the interval -1 < s < 0. Equation (1) is nonlinear and of the 



Volterra type (ref. 2). If taken as an equation for r(s) , it has a singularity at s = 0. 
Equation (2) is linear and of the Fredholm type. It hac u singularity at q = s. 

The solution of the problem depends on the proper handling of these two singulari- 
ties. The given conditions on r(s) and d(s) are 

t(s) = t(-s) with t(1) = 0 

6( s) = -6(- s) with 0(0) = 0 

where r(s) is an even function of s, and 8(s) is an odd function of s. 

METHOD OF SOLUTION 

The numerical solution is obtained by using Picard's method of iteration. Equations 
(1) and (2) are rewritten as 



The next estimate 0 K is defined by 

0 K (s) = e K _ 1 (s) + jSD(s) (5a) 

where D is the suggested correction to defined by 

D(s) = U K (s) - 0 K _ 1 (s) (5b) 

and £ is the relaxation parameter. 

The three steps of the iterative loop are briefly 
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(1) Solve equation (3) for t^-(s) . 

(2) Solve equation (4) for U^-(s) . 

(3) Combine 0^ ^ and U K by equations (5a) and (5b) to get a new estimate 0^. 

The starting estimate used was 

*o< s )=— (6) 

0 4 

Two aspects of the iterative process will be discussed in greater detail in the follow- 
ing sections: 

(1) In the section INTEGRATION TECHNIQUE, the principal value 

(2) In the section ACHIEVEMENT OF CONVERGENCE, convergence of the iterative 
process by 

(a) Proper handling of the singularity in equation (3) 

(b) Proper choice of the relaxation parameter /3 of equation (5a) 


INTEGRATION TECHNIQUE 


To evaluate U K from equation (4) requires a technique for dealing with the princi- 
pal value of an integral. Since the singularity is at q = s, and t k is known at all s 
points, -1, -1 + As, . . . , 1 - As, and 1, the integral is split up into at most three parts, 
with the second, containing the singularity: 

U K = I 1 + I 2 + I 3 


•S-As r S+AS fl. 00 

+ J + J 

-1 s-As ^s+As 


The 1^ and 1^ integrals may be evaluated by Simpson's rule. The Ig integral is 
treated as follows: 

For s = 0, As, 2As, . . . , and 1 - As expand r^- in a Taylor series around q = s 
and truncate after three terms to get the approximation r* : 

T* (q) = t k (s) + Tg-(s) (q - s) + T^(s) s t ~ (8) 

2 

Replace equation (8) by its finite difference analog to get the approximation r n (q): 

a 
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r a (q) - t k (s) + 


r K (s + As) - t k (s - As) 
2As 


(q - s) 


t k (s + As) - 2r K (s) + t k ( s - As) ^ _ g ^2 


(9) 


As" 


Replace r K (q) in I 2 by r R (q) to get 


I 2 (s) = 


t k (s) 


/iS+As 

^s-As 


An T t^(s + As) - T v (s - As) /-S+As 

M- + * / dq 

q - s 2AS7T ‘'s-As 


t k (s + As) - 2r Tr (s) + r^-Cs - As) 


K v ' T 'K v 


7tAs 


J aS+As 

s-As 


SLLSdq 

2 


( 10 ) 


The principal value of the first integral is zero. The second integral is 2As, and 
the third integral vanishes to yield 


Tr^ (s + As) - Tt^-(s - As) 

I 2 (s) = — for s = 0, As, . . . , 1 - As 


(ID 


At s = 1, I 2 reduces to 


f 

‘'l-As 


r K (q)dq 

q - 1 


For this point, all derivatives appearing in the Taylor series expansion of T K (q) 
about the point q = l are one-sided. Use T. to denote 7 K (1 - jAs); then, since 
T q = 0, the equations for the finite difference analogs of r’ and r" become 


rfed) = — +°.5 

^ As 



( 12 ) 
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and 


~2T 1 + T 9 /-2T 1 + T 9 T 1 - 2T„ + T„ 
r^(l) = 1 + 12 1 2 3 


As 


As 


As 


■ 5T 1 + 4T 2 - T 3 
As 2 


(13) 


Substitute equations (12) and (13) in equation (8) for s = 1 to get the approximation 

T a (q) ; 


T - 4T. -5T + 4T 9 - T, ( 

r Q (q) = -2 — (q - 1) + * 2 - (q_nl)_ 


2As 


(14) 


As^ 


This formula is used in the equation for 


i 2 (« = 


1 f T a (■»>*» 

7T / q - 1 
“l -As 


(15a) 


to yield 


yi) 


-3T t - 2T 2 + Tg 
An 


-3r^(l - As) - 2 t k -(1 - 2As) + 7>(1 - 3As) 

= - - (15b) 

An 

The evaluation of 1^ and Ig was handled first by Simpson’s rule formulas. Later 
a second method was applied. Both methods yielded sufficient accuracy for this problem 
and for the test case. But other problems involving principal value integration may bene- 
fit from the more careful technique of the second method. Therefore, its derivation is 
given here . 

Note that the 1^ integrand r K (q)/(q - s) is large and of one sign near s - As, 
whereas the I 3 integrand is large and of opposite sign near s + As. Using differences 
of large numbers may increase the truncation error. Combining symmetric portions of 
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L and I„ should reduce the truncation error. 
1 3 

Let S. be the integral sum 


S. = 
] 


/ s-jAs 
(j+l)As 


T K (q)dq 

q - s 


/ s+(j+l)As 
oi-jAs 


T K (q)dq 
q - s 


where j = 1, 2, 3, . . . , j max and 


^max ^ 


( 16 ) 


( 17 ) 


where i is such that s = (i - l)As and N is the number of points in O^ssl such 
that As = 1/(N - 1). Then the sum of 1^ and Ig becomes 


I x (s) + I 3 (s) = 


I 


2s-l 


T K (q)dq 

q - s 


max 


S, 


j=l 


(18) 


where much of the important cancellation of large numbers occurs in the computing of 
S. for low j. 

The S. are computed as follows. Make use of a change of variables to obtain for 
equation (16) 


S. = 
] 




f 2 (x)dx 

x 


(19) 


where a ^ < a 2 . 

Fit a parabola through the three values of f^ at x = ~a 2 , 

f^(x) = a^ + b^x + c^x 2 

and fit a parabola through the three values of f 2 at x = 2a ^ - 

2 

f 2 (x) = a 2 + b 2 x + c 2 x 


-a-p and o? 2 - 2 o'., to get 

( 20 ) 

o; 2 , a? ^ , and o! 2 to get 

( 21 ) 
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Then S. 
] 


may be integrated in closed form to get 


S j = a i la 


ot i 


a< 


+ b l (a 2 “ «1> + V 


a. 


or 


a 2 

+ a„ In — + boCoJo 
“1 




( 22 ) 


s j - (a s 


0/n (Cn ~ Cj 

a^) In — + + b 2 )(c* 2 - O'j) + — — 

0 !^ 2 




(23) 


As a test, this integration process was applied to the function 



(1 - q 2 )dq 
q - s 


(24) 


The numerator of the integrand in equation (24) is like r in being even and zero at 
q = 1. The exact value of g(s) is 


g(s) = (1 - s 2 )ln — - - 2s 
1 + s 


(25) 


Note that g(l) = -2. 

The following table gives the integral of equation (24) for values of s together with 
the correct value of g(s). The integration technique is clearly accurate. 
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Distance 
parameter , 
s 

Integral of eq. (24) 

Numerical value 

Correct value 

0.1 

-0.39866259 

-0.39866398 

.2 

-.78924440 

-.78924650 

.3 

-1.1633240 

-1.1633257 

.4 

-1.5117280 

-1.5117302 

. 5 

-1.8239566 

-1.8239592 

.6 

-2.0872253 

-2.0872283 

.7 

-2.2846440 

-2.2846465 

.8 

-2.3909990 

-2.3910008 

.9 

-2.3594423 

-2.3594434 

.99 

-2.0853368 

-2.0853368 


ACHIEVEMENT OF CONVERGENCE 


One action necessary to achieve convergence is to handle properly the singularity in 
the equation for r at s = 0. In equation (3), the term in brackets can be written as 


x-VT7 




-3 

Its cube behaves like s near s = 0. 

Since the upper integral on the right side of equation (3) behaves like 


(26) 


s rs 3 

q sin(Kq)dq ~ / Kq 2 dq = (27) 

Jo 3 

o o 

it becomes obvious that the right side of equation (3) behaves like s /s in the vicinity 
of zero. 

Because of the singularity at zero, the numeric difficulties become more severe as s 
approaches zero. Any errors generated in the evaluation of the upper integral are mag- 

3 

nified when divided by s and lead to instability. Since t is known to be an even func- 
tion of s with zero derivative at s = 0, it seems logical to impose proper behavior in 
the neighborhood of zero numerically. 

Another necessary action to achieve convergence was to choose the proper relaxation 
parameter /3. (See eq. (5a) for its definition.) 
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The use of j3 = 1.0 as the coefficient of the suggested change led to an oscillation 
of about 5 percent in the value of t( 0) (see fig. 1). Figure 2 shows how the oscillation 
drops to less than 3 percent when /3 is 0. 9. Figure 3 shows how the oscillation drops 
to less than 0. 5 percent when /3 is 0. 7. Figure 4 shows that with /3 = 0. 5 the oscilla- 
tion disappears, and convergence is very good after 18 iterations. 

The converged shape of the r(s) function is shown in figure 5. The converged shape 
of 6{s) is shown in figure 6 . 


CONCLUDING REMARKS 

This report shows that there were two necessary steps required to achieve the solu- 
tion of the original system of two integral equations . The steps taken were 

(1) The T-g(s) function was isolated in the neighborhood of s = 0, and a smoothing 
process was applied to remove the erratic behavior of r^- near s = 0. 

(2) A relaxation factor equal to about 0. 5 was used to remove the oscillations of the 
r K^(s) curve. 

This method, then, provided an accurate solution to the system of two integral equa- 
tions, as is shown in the section of the report on accuracy. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 10, 1971, 

129-04. 
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Figure 4. - Behavior of t(0) 
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